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ABSTRACT 

We study the evolution of a magnetic arcade that is anchored to an accretion disk and 
is sheared by the differential rotation of a Keplerian disk. By including an extremely low 
f^ ' external plasma pressure at large distances, we obtain a sequence of axisymmetric magnetostatic 

C^ , equilibria and show that there is a fundamental difference between field lines that are affected by 

the plasma pressure and those are not (i.e., force- free). Force- free fields, while being twisted by 
the differential rotation of the disk, expand outward at an angle of ~ 60° away from the rotation 
axis, consistent with the previous studies. These force-free field lines, however, are enclosed by 
(^ ' the outer field lines which originate from small disk radii and come back to the disk at large 

CN \ radii. These outer fields experience most of the twist, and they are also affected most by the 

external plasma pressure. At large cylindrical radial distances, magnetic pressure and plasma 
pressure are comparable so that any further radial expansion of magnetic fields is prevented 
^S| ' or slowed down greatly by this pressure. This hindrance to cylindrical radial expansion causes 

most of the added twist to be distributed on the ascending portion of the field lines, close to 
j^ I the rotation axis. Since these field lines are twisted most, the increasing ratio of the toroidal 

Bff) component to the poloidal component -B_r,z eventually results in the collimation of magnetic 
energy and flux around the rotation axis. We discuss the difficulty with adding a large number 
of twists within the limitations of the magnetostatic approximation. 

i-G ■ Subject headings: Accretion, Accretion Disks — Magnetic Fields ~ MHD — Plasmas 

O I 1. INTRODUCTION progress has been made recently in the hydromag- 

^ : The process of forming colhmated jets/outflows gf nTrlTrP^'^VTosIr' D^n'^mi^ MHo'tJ"^ 

C^ ■ due to disk accretion onto central compact objects Yl ,, 

■ ■ is thought to depend on how magnetic fields be- 1^*^°"^ ^^ ^^^ near jet region have been carried out 
r1 ^ have when they are swirled around by the accre- ^^ severa groups (Beh 1994, Ustyugova et al. 1995, 

^ ' ,- I- , rpi t J 4- J- 4-u- 1999; Koldoba et al. 1995; Romanova et al. 1997, 

^Vj , tion disk. The progress of understanding this pro- '^ , „^ ,ot-.t ^, 

d ■ cess has, however, been hindered by the significant l^^^'^ ^^f «* al. 1997; Ouyed & Pudritz 1997; 

■ ■ ' lack of knowledge on the global magnetic field con- Krasnopolsky, Li, & Blandford 1999) The Simula- 

figuration in/near the accretion disk (see Okamoto *^°^ ^^^^^ °f Ustyugova et al. (1999 m the hydro- 

1999 for detailed critiques on many MHD models; see magnetic wind regime indicates that the outflows are 

also Blandford 2000 for a recent review). An ordered accelerated to super-fast magnetosomc, super-escape 

magnetic field is widely thought to have an essential speeds close to then region of origin (-80 x the m- 

role in jet formation from a rotating accretion disk. "^^^ ^^^^^s o^ ^^^ disk), whereas the collimation oc- 

Two main regimes have been considered in theoreti- ™^^ ^^ much larger distances. Poyntmg flux mod- 

cal models (see Lovelace, Ustyugova & Koldoba 1999 ^^^ ^^^ *^^, origin of jets were proposed by Lovelace 

for a review), the hydromagnetic regime where the en- (^9^6) and Blandford (1976) and developed further 

ergy and angular momentum is carried by both the ^y^^,^^^^^' ^^^S, and Sulkanen 1987 Lynden- 

electromagnetic field and the kinetic flux of matter, Beh (1996), and Colgate & Li (1999). In these mod- 

1 ,1 D ,■ zj ■ 1, j-u A els the rotation of a Keplerian accretion disk twists 

and the Foynhng flux regime where the energy and , ■ , i r. i i i i- i i- t i i ■ 

angular momentum outflow from the disk is carried -a r v,^-it,-v,' 

predominantly by the electromagnetic field. Major j & 
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mentum (in the twist of the field) and energy (in the 
Poynting flux) away from the disk, thereby facilitat- 
ing the accretion of matter. Most recent computer 
simulations using the full axisymmetric MHD equa- 
tions have been in the hydromagnetic regime. How- 
ever, recent simulation studies have found jet outflows 
in the Poynting flux regime (Romanova et al. 1998; 
Ustyugova et al. 2000). 

One class of models deals with a simplified limit, 
where magnetostatic and force-free conditions are as- 
sumed. This kind of problem is very clearly stated in 
the introduction of Lynden-Bell & Boily (1994, here- 
after LB94; see also Lynden-Bell 1996): A field rooted 
in a heavy conductor on z = (i.e., the rotating disk) 
pervades a perfectly conducting force-free medium in 
the region z > 0. On the disk surface {z = 0), the 
field passes upwards from the region < R < Rq and 
returns downwards in the region Ro < R < -Rmax- 
The disk is now rotated about its axis according to 
the Keplerian motion, Q{R) oc R~^''^. The problem 
is to determine the magnetostatic field configuration 
when the disk has gone through a specific number of 
turns. They found that fields, instead of collimat- 
ing along the rotation axis, expand along an angle of 
K, 60° away from the rotation axis. 

The behavior of twisted magnetic field lines that 
are anchored in a perfectly conducting medium has 
originally been considered in the solar corona context 
(see, e.g., Aly 1984, 1991, 1995; Mikic et al. 1988; 
Finn & Chen 1990; Sturrock et al. 1995). The di- 
rect application of those studies to accretion disks 
has been fairly recent (e.g., Appl & Camenzind 1993; 
Konigl & Ruden 1993; LB94; Lynden-Bell 1996; Bar- 
dou & Heyvaerts 1996; Goodson et al. 1999; Uzden- 
sky et al. 2000). The important role of an external 
plasma pressure was first emphasized by Lynden-Bell 
(1996) where he argued that field lines can not sim- 
ply expand to infinity as shown by LB94 because the 
required PdV work becomes too large. He then con- 
structed a simplified cylindrical model where a mag- 
netic helix/jet is bounded at large radius by exter- 
nal plasma pressure. But the question left open is 
whether the plasma pressure can indeed play a con- 
fining role in a self-consistent global treatment of the 
twisting and expansion of magnetic fields driven by 
the disk shear rotation. 

In this study, we obtain self-consistent global so- 
lutions of axisymmetric magnetostatic configurations 
by strictly enforcing the Keplerian shear condition on 
the flux surfaces that are anchored in the disk. The 
assumptions we use are explained in detail in §g. We 
formulate our problem in §^, along with the relevant 
equations and methods to solve them. Results are 
given §^ and final conclusions and discussions are in 



2. BASIC ASSUMPTIONS 

The equation of motion in the non-relativistic ideal 
magnetohydrodynamics (MHD) limit is simply 



p-rr = - J X B + pg 
at c 



VP 



(1) 



where p is the plasma density, v the flow velocity, J 
the current density, g the gravitational acceleration, 
and P is the plasma pressure. 

We restrict our attention to axisymmetric magne- 
tostatic configurations. The underlying assumption 
is that the speed of the field line twisting by the disk 
is slow so that the system quickly reaches an equilib- 
rium. We can then treat the evolution of field con- 
figurations as sequences of magnetostatic equilibria. 
We make the further assumption that the magnetic 
field plays the dominant role with plasma flow v and 
gravity pg much less important. Then, the steady 
state equation is 



J X B = c VP 



(2) 



The main astrophysical motivation for keeping this 
pressure term is that at sufficiently large distances 
away from the disk, the plasma pressure will be com- 
parable to the magnetic pressure, thus becoming dy- 
namically important. When VP — > 0, we have the 
so-called force- free (FF) limit. 

We now discuss the physical conditions where these 
assumptions apply. There are at least two relevant 
speeds in this problem which are related to two phys- 
ical aspects. The first one is the velocity of Keplerian 
disk rotation, uk, which describes the rate of field line 
footpoint movement. The second is the relaxation 
speed, vx, for magnetic fields to reach equilibrium, 
which is usually achieved by MHD waves going back 
and forth in the system. This is essentially the Alfven 
speed VA and is ultimately limited by the speed of 
light c. Together with these speeds, there are two 
relevant timescales: the disk rotation period Tk = 
27rPmm/^K5 and the relaxation time Tx = imax/^Xj 
where Pmin is the inner disk radius, and L^ax the sys- 
tem dimension, respectively. In order to ensure the 
steady state condition on the timescale of Tk, one 
requires that L^ax <■ {vx/vK)Rinm- 

In order to simulate a system that is much larger 
than Pmiri) these considerations indicate that in a 
black hole accretion disk system, field lines close to 
the black holes are likely twisted too rapidly to be 
treated by the magnetostatic equations. In other 
words, the dynamic pressure from the inertial term 
pdv/dt in equation (||) is likely to be very important 
during the expansion, thus the conclusions from this 
study might not apply in this limit. On the other 



hand, for field lines further away from the black hole 
and/or for accretion disks around systems like young 
stars, the rotation speed is small enough that the 
magnetostatic limit could still apply. 

Another relevant process is the footpoint drift 
across the disk. Magnetic field threading the disk 
tends to be advected inward with the accretion flow, 
but at the same time it may diffuse through the disk 
owing to a finite magnetic diffusivity rim of the disk. 
The outward drift of the magnetic field in the disk 
occurs at speed Ur = {rjm / h) tan{9) , where h is the 
half-thickness of the disk, tan(0) = (^Br/Bz)z=Oi and 
a smaller, second order diffusion term {'qmd'^Bz/dr'^) 
has been omitted (Lovelace et al. 1997). For cases 
where the diffusivity is of the order of the viscosity 
and where the viscosity is given by the Shakura and 
Sunyaev (1973) prescription u = acgh (with a < 1 
and Cs the midplane sound speed), the diffusive drift 
speed is Ur ~ aCgianiO). This is larger than the ra- 
dial accretion speed Vr ~ —aCs{h/R) due to viscosity 
alone. But what is important here is that the field 
drift and the accretion speeds are much less than the 
Keplerian velocity of the disk for Cs <C Vk- For this 
reason the fields can be treated as frozen into the disk. 
This point was made by Ustyugova et al. (2000) and 
also discussed by Uzdensky et al. (2000). 

3. BASIC EQUATIONS 

We have formulated the problem discussed in ffl in 
both cylindrical {R, (p, z) and spherical (r, 6, (p) coor- 
dinates. Axisymmetry is assumed in both cases. The 
magnetic field can be written as: 



and for spherical coordinates, 



B = Bp + 



V* X V(^-h 



(3) 



where the poloidal component Bp is in the {R, z} or 
{r, 9} plane, 2tt^ is the total poloidal flux through 
the disk, and contours of ^ label the poloidal field 
lines. The poloidal magnetic fields are 

R - ^ ^"^ R - ^ ^* (A\ 



for cylindrical coordinates and 

1 a* 



Br 



Bn 



1 a^ 



r^ sin 9 do ^ " r sin 9 dr 

for spherical coordinates. The current density 

—J = -A*^V0 + V{rB^) X V(j) , 



(5) 



(6) 



where the operator A* is, for cylindrical coordinates, 



A*^ = 1 

dR'^ RdR dz^ 



(7) 



At^ 



d^ 



+ 



1 ^2 



1 



d 



Qj.2 J.2 QQ2 J.2 ^qt^ Q QQ 

Equation (^) implies that 

B-VP = =^ p = p(^^). 



(8) 



(9) 



In other words, the gas pressure is constant along field 
lines. In cylindrical coordinates, equation (||) can be 
written as the well-known Grad-Shafranov equation, 

dP 
A*^V^ + RB^V{RB^) = -AttR^—V^ , (10) 

from which we define RB^j, = H{'^) and we get 

A^* + d[H^/2]/d'^ + AirR^dP/d-^ = . (11) 

Similarly, in spherical coordinates, we define 
r sin 9B^ = H{^) and 

A^^ + d[H'^/2]/d^ + 47r(r sin OfdP/d^ = . (12) 

In our spherical calculations, we actually use ^ = In r 
to concentrate the uniform grid in ^ for better reso- 
lution at small r. Equation ( p!^ ) then becomes 

A^^^ + r'^d[H'^/2]/d^ + 4.i:{r'^ sm9'^)dP/d^ = , 

(13) 
where 



a:,* 






d_ d^ 
d^ ^ 09^ 



1 d 

tan 9 89 



(14) 



The quantity H{^) = RB^{R, z) is (2/c) times the 
total current flowing through a circular area of radius 
R (with normal z) labeled by ^{R,z)= const. The 
development of the toroidal field component from an 
initial purely poloidal field comes from the differen- 
tial rotation of the disk onto which footpoints of the 
same field lines are anchored. 

3.1. Boundary Conditions 

Since there is no direct observational information 
on how magnetic fields are distributed on the sur- 
face of an accretion disk, we have made the follow- 
ing assumptions. If the magnetic fields on the disk 
is initially generated via an accretion disk dynamo 
that makes a quadrupole or dipole-like field, magnetic 
fields will emerge out of the disk at smaller radii and 
go back to the disk at large radii. We use Ro to repre- 
sent the 0-point of the field in the disk, i.e., the radius 
where i?^ = on the disk. Furthermore, in a fully 
self-consistent treatment including the back-reaction 



on the disk, Rq also marks a separation location in- 
side which angular momentum is lost and transmitted 
to the outer part of the disk by field line tension. 

We study the expansion of this cylindrical mag- 
netic arcade anchored on the disk. We further as- 
sume that the overall strength of \Bz\ decreases as 
a function of radius, except near the 0-point where 
\Bz\ = 0. This is roughly consistent with the fact 
that the thermal pressure of the disk (which anchors 
the fields) decreases radially as well. Specifically, we 
have chosen a computational domain that is a box 
with < -R < -Rmax and < z < z^^ax in the cylin- 
drical case and/or a sphere with rmin < r < rmax in 
the spherical case. The outer boundary is assumed 
to be perfectly conducting (^ = 0). At z = (along 
the disk), we assume that the disk is perfectly con- 
ducting as well. An initial poloidal field distribution 
is specified ^(i?,0) [or *(r,7r/2)] as 



fl2 

*(i?,0) oc { R'^ 



ii° 



for R < Re ; 

for Rc<R< Ro] (15) 

for Ro < R< Rout, 



where we have joined these parts smoothly. Fig- 
ure |l| shows the distribution of ^ and \Bz\ on the 
disk surface. The poloidal flux is distributed be- 
tween -Rmin ~ 10~^ and Rout ~ 0.02. A core radius 
Re ~ lOfimin is uscd, insidc of which the B^ compo- 
nent approaches a constant. The index a is chosen 
to be 5/4, so that B^ is decreasing as i?"~^, revers- 
ing its sign at the 0-point Ro ~ 0.01, and continuing 
decreasing until -Rout- Both ^(R,0) and \Bz\ remain 
zero between i^out and -Rmax = 1- These choices for 
Ro and Rout ensure that the initial magnetic field is 
distributed far inside the outer boundary i?max- We 
have tried many different initial field configurations 
and found that our main conclusions do not depend 
on these particular choices of parameters, as long as 
magnetic fields stay away from the outer boundary 
-Rmax or Tinax- The dependence on ^{R,z = 0) near 
the 0-point is very weak since those field lines are 
hardly twisted at all. 

3.2. External Plasma Pressure 

The physical picture we have in mind for account- 
ing for the pressure of an external plasma is that mag- 
netic field close to the disk is force-free with little in- 
fluence from the plasma pressure. At large distances, 
however, the magnetic pressure decreases sufficiently 
so that plasma pressure becomes important. In other 
words, the magnetic fields are bounded above and on 
the sides by an ambient medium which is perfectly 
conducting and mostly devoid of any magnetic flux. 
This forms a conducting, in general, "deformed box" 
inside of which the magnetic field dominates and the 
gas pressure is very small. 



To model the external plasma pressure, we adopt 
the dependence 

P(*)=Peexp[-(*/^,)2] , (;Lg) 

so that 



dP/d^ = -{2PJ^I) ^ exp[-(*/*c)^] 



(17) 



where Pc and ^c are the input parameters. For most 
of the results presented here, we have chosen Pc = 0.1 
and ^c = 10~^. The physical meaning of Pc is that its 
magnitude gives an estimate of the overall importance 
of plasma pressure. Using the example given in equa- 
tion (p!5|), the maximum Bz near Rmin is ~ 6 x 10^, 
(see Figure ||) , so the ratio of the maximum magnetic 
pressure to the maximum gas pressure in the compu- 
tational domain is ~ 10^^. The quantity ^c gives the 
fraction of poloidal fiux that is affected by the plasma 
pressure. 

The prescribed plasma pressure (equation |l^) is 
kept fixed at all times throughout the sequence of 
equilibria. This is because the prescribed pressure 
here is meant to mimic a constant external pressure 
boundary (such as pressure from interstellar or in- 
tergalactic medium) that is reacting to the "push" 
by magnetic fields. This interpretation is precise in 
the limit of ^c — ^ 0. We have used a very small 
but finite ^c for the convenience of modeling both 
the Lorentz force and pressure gradient terms simul- 
taneously, rather than having to specify the pressure 
gradient term through a boundary condition, espe- 
cially when the shape of this boundary is determined 
dynamically by the push from the magnetic fields and 
is not known a priori. The pressure in the flux tubes 
with ^ > ^c, however, might change due to the flux 
volume expansion. If the entropy of each flux tube is 
conserved, then pressure has to decrease with expand- 
ing volume. On the other hand, if there is enough 
heat flux between the disk and the flux tubes, then 
the entropy of a flux tube is not conserved (Finn & 
Chen 1990). Either way, since the pressure in the flux 
tubes with ^ > ^c is already set to be exponentially 
small (equation ^) initially, any further decrease in 
pressure will not alter our conclusions. 

We make two more observations about the plasma 
pressure term. First, equation ( [I^ ) gives that, when 
^ — > (i.e., at the boundary), the pressure P — > Pc 
and its gradient dP/d^ -^ 0. Both features are phys- 
ically plausible. Second, it is interesting to note that 
the pressure effect enters equations ([ll|) and (|l3|) with 
a geometric factor R^ or (rsin^)^. This factor actu- 
ally comes from the J x B term. It is only under 
the equilibrium condition can it be "moved" to be 
a multiplier of the pressure term. Consequently, the 
pressure tends to prevent expansion away from the 
rotation axis but has relatively little constraining ef- 
fect along the rotation axis. 



3.3. Input Keplerian Field Line Twist 

From the distribution of ^{R, 0) and the Keplerian 
rotation ^{R) oc R^^''^, we can define the required 
twist on each field line A<^(^). The twist of a given 
field line going from an inner footpoint at Ri to an 
outer footpoint at R2 is proportional to the differen- 
tial rotation of the disk. Using a cylindrical coordi- 
nate, for a given field line, we have 



Rd4)/B^ = dip/B, 



P ! 



where dtp = s/dR^ + dz"^ is the poloidal arc length 
along the field line, and Bp = \l B% + B^. Then the 



'R 



total twist of a field line is 

/■2 B, 

I dip- 



A^fl-) 



30 



H{^) 



RBp "' ' ji R'^Bp 



(18) 



We denote the quantity V'{^) = J^ dip/R^Bp. The 
field line twist after a time t is 



A$(^) 



VLo t 



Ro\'/' 

rJ 

(fto t) H"^) 



Ro\'/' 
R2 



(19) 



where Qq = \/GWJR^ is the Keplerian angular fre- 
quency at Ro of an object of mass M, and J-' is di- 
mensionless. The rotation profile is assumed to be- 
gin deviating from Keplerian and approach a con- 
stant when R < Re (the bottom panel of Figure 0). 
There exists a region, however, between ^ = and 
^min = "^{Rmin) where the twist decreases as ^ de- 
creases (i.e., getting closer to the z-axis and the sur- 
rounding wall). These field lines are not explicitly 
followed in our calculations. But since -ff(^) — > 
as \I' — > 0, the twist A$(^) approaches zero as well 
[equation ([l8|)]. 

3.4. Method of Solution 

We solve equations (0) and ([T3|) following the ap- 
proach given in Finn & Chen (1990), which employs 
several levels of iterations. For the very first step, 
we will use a trial H(^), then we solve equation 
(PI ) or (|l^) using the Successive Over-Relaxation 
(SOR) method (see for example. Potter 1973) for 
^^{r,9) -^ ^''~^^{r,9). From ^''+^(r, 61), we integrate 
along the field lines to obtain V'{'^). We typically 
trace ~ 200 field hues with 3 x 10"^ < ^ < 1 for 
calculating the poloidal current profile H{^). Us- 
ing equation (18), we update H('^) using the input 
A<l>(^). This procedure is then repeated. The ad- 
vantage of the simple SOR solver is that the nonlin- 
ear source terms d{H'^)/d^ and dP/d'if can be easily 
included in the iterations so that convergence can be 



achieved fairly quickly (with a typical relative residue 
less than 10"^) 

In summary, we compute global magnetostatic so- 
lutions of ^(r, 9) using the Grad-Shafranov equation 
in axisymmetry by requiring that the twist on each 
field line follows a specified function A$(^), derived 
based on the Keplerian disk rotation. 

4. RESULTS 

We have made many runs with different choices of 
the initial field configuration on the disk ^{R, 0), the 
twist profile A<1>(^), and the pressure profile P{^)- 
Runs are made in both cylindrical and spherical coor- 
dinates. We present most of our results using the fol- 
lowing parameters: A In r— spherical coordinate sys- 



tem is used, with r^ 



10" 



1. The radial 



and 0— angle grids are 367 x 51. The initial field con- 
figuration is the same as shown in Figure |^, where 
Ro = 0.01 and Rout = 0.02. The smallest ^ value of 
the field lines we track is ^min = ^{Rmin) = 3 x 10~^. 
The plasma pressure parameters are Pc = 0.1 and 
^c = 0.001. The maximum \Bz\ on the disk is 
~ 6 X 10^, which implies a ratio of ~ 10^^ between 
the maximum magnetic pressure and the plasma pres- 
sure. (Even larger pressure ratios can be achieved 
with relative ease.) To indicate the progressively in- 
creasing twist added to the field lines, we use the 
notation time "i" [equation (|l9[)] in units of revolu- 
tions of the inner most flux line (^min) around the 
z-axis, i.e., A^^ax = A$(^min) = 27rt. 

Figure |2| describes our overall physical picture with 
two different physical regimes. For < \1' < ^o the 
plasma pressure is dominant. This region is labeled 
as "plasma-pressure" (PP). For ^c < ^ < 1; the 
magnetic pressure dominates, the region is force-free 
(labeled as "FF"). So, we have effectively made two 
ideal MHD fluids, one (PP) has a plasma /3 — > 00 and 
the other (FF) has /? = 0. The key question we are 
addressing is: "how does the shape of ^ = ^c bound- 
ary change while field lines are being twisted by the 
Keplerian disk rotation?" In other words, the bound- 
ary between the PP and FF regions evolves according 
to both the expansion and "pushing" from the FF re- 
gion and the "hindrance" of the plasma pressure. 

Although some details may differ, we generally find 
that magnetic fields evolve in the following manner: 
(1) in the FF region, field lines expand primarily to- 
wards the large radius along an angle (from the z- 
axis) of (9 ~ 45° - 60°; (2) in the PP region, field 
lines expand predominantly vertically (along the z- 
axis) with some slight radial expansion; (3) more 
twist causes further expansion along the z-axis but 
subject to a break-down of the magnetostatic, equi- 
librium assumption. We now discuss these features 
in detail. 
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4.1. Equilibrium Sequence with Increasing Twist 

Figure |3| shows the evolution of magnetic fields in 
the poloidal plane as twist is added, with t = 0, 1, 2 
and 4 turns, respectively. The field lines shown are 
evenly spaced in log;^o(^)) with the outer- most line 
having *& ~ 10~^ and the inner-most line having 
^ ~ 0.5. Here, we use the term "outer" to refer to 
field lines in the PP region, which originate from the 
smallest radii and return back to disk at the largest 
radii. The term "inner" refers to field lines around 
the 0-point in the FF region. Consequently, because 
the field lines are twisted according to the Keplerian 
rotation (largest rotation at the smallest radius), the 
outer field lines show substantial "movement" due to 
the added twist, whereas field lines around the O- 
point show little change. 

The corresponding results for the poloidal current 
H{^ and the resultant twist A<1> are plotted in Fig- 
ure ^ It is interesting to note that -ff (^) can be ap- 
proximated by two power-laws of ^, where H{^) oc ^ 
when ^ <^c and flattens somewhat for ^ > ^c- The 
twist shown in the lower panel is derived by integrat- 
ing along each field line. They indeed follow closely 
the input profile due to Keplerian rotation. 



4.2. Role of External Plasma Pressure in 
Collimation 

The striking feature of the sequence of equilibria 
in Figure y is the way the outer field lines expand. 
To best understand this behavior, we separate the 
field lines into two groups according to the ratio of 
^/^c- These two groups have fundamentally differ- 
ent expansion behavior, which is illustrated in Figure 
III, where the expansion of three different field lines 
(^ = 10^2 10~3,10~'^) is shown for f = (solid), 
t = 2 (dotted), and i = 4 (dashed), respectively. 

Field line expansion in the FF region has been stud- 
ied in detail in many previous studies (see Introduc- 
tion). The Lorentz force is the only force available 
(i.e., the magnetic pressure gradient is balanced by 
the tension force). As found for example by LB94, 
the field lines expand to large radii along an angle 
of ~ 60° from the z— axis. This is consistent with 
the left panel of Figure ^ (also the inner field lines in 
Figure ^) . 

The field line expansion in the PP region, however, 
is fundamentally different from the FF regime. This 
is because the magnetic pressure becomes very small 
at large radii and eventually at some distance it be- 
comes comparable to the surrounding plasma pres- 
sure. At this location, further expansion of the mag- 
netic fields is greatly hindered by the plasma pres- 
sure. At the same time, increasing the twist acts 
to increase the ratio of B^/Br. Consequently, these 



field lines are increasingly "pinched" around the z- 
axis, and they eventually expand along the z-axis in 
a "collimated" fashion, as shown by the middle and 
right panels of Figure |5| and the late time in Figure 

The middle panel of Figure |5| is particularly impor- 
tant. It shows how the critical boundary at ^ = ^c 
(see Figure ^) evolves with increasing twist. Note 
that its initial shape is quasi-spherical, but a clear 
collimation along the z— axis has developed by i = 4. 

A small amount of poloidal flux (with ^ < ^c) is 
in the PP region and is twisted as well. One can ask 
whether these outer field lines have contributed im- 
portantly to forming the new shape of the ^ = ^c 
boundary. To answer this question, we have per- 
formed runs with a very small A<1>(^) for ^ < ^cj 
instead of the near constant A<1>(^) presented in Fig- 
ure ^. We find that as long as A$(^c) remains the 
same, its shape does not depend on the twist profile 
for ^ < ^c- This means that the small amount of 
magnetic flux beyond the ^ = ^c boundary does not 
change our results noticeably. 

We have also performed runs with various magni- 
tudes of the plasma pressure. For the purely force- 
free case, we confirm the results by LB94 that the field 
lines expand to large radii along an angle of ~ 60° 
from the z— axis and only a small total twist (less than 
half a turn) can be added before the fields are packed 
against the outer boundary. With a finite pressure, 
however, the collimation effect is always observed as 
long as the simulation region is large enough or the 
pressure is large enough so that magnetic fields are 
not directly packed against the outer boundary. 

Thus, the external plasma pressure, no matter how 
small it may be, plays a fundamental role in causing 
the field lines to collimate along the rotation axis. It 
essentially stops (or greatly slows down) the radial 
expansion of the field lines, and in the meantime al- 
lows the build-up of the toroidal component. This 
interpretation is also consistent with the fact that 
field lines are closely packed at large radius (stopped 
by plasma pressure) but are loosely spaced along the 
z-axis (see the "i = 4" panel in Figure 3|). 

The results shown in Figures |3| and |5 indicate that 
the solution suggested by Lynden-Bell (1996) is in- 
deed possible (at least in the early expansion stage 
of the helix formation), even though details are dif- 
ferent. Here, we have shown global self-consistent 
solutions where the field lines are twisted according 
to the disk Keplerian rotation. 

4.3. Distribution of Twist and Magnetic Energy 

Figure ^ shows how its twist (A0) is distributed 
along this field line as a function of z at t = 4. 
We have plotted another field line with a larger 



^2 = 10^^ at t = 4 as well. There is an impor- 
tant difference between these two field lines in the 
distribution of the twist A<I>(^). For ^i, most of its 
twist (> 75%) is distributed on the ascending portion 
of the field line and near the z-axis; the rest is dis- 
tributed while the field line spirals down back to the 
disk. For ^2) however, the twist is distributed uni- 
formly between ascending and descending portions. 
To understand this difference, we can write 



Hi^) 



d(t) 
dz 


= 




= 



Rd^/dR 



Note that H{'^) is constant along a field line. Imagine 
a small flux tube (^ ^ ^ + d'^) originating from the 
inner region of the disk. Since the field line remains 
at small R for the ascending section in the presence of 
plasma pressure, B^R? ~ 2^. This is the reason why 
d(f)/dz following a flux tube is roughly constant for the 
ascending part of the flux tube. For the ^i flux tube, 
however, its expansion at large radius (the descending 
portion) is strongly constrained by the plasma pres- 
sure, so that the cross section area for the descending 
fiux tube is much smaller than it would have been 
without the plasma pressure. Consequently, its B^ 
component on the descending portion of the 'J' = ^i 
surface varies more slowly than i?~^, so that the rate 
\d(f)/dz\ is smaller. Overall, the twist is then non- 
uniformly distributed on a fiux line. (See Ch .9 of 
Parker 1979). Thus, the presence of plasma pressure 
leads to a strong collimation and a concentration of 
the twist to the collimated region. 

The total injected toroidal flux !F^ by the disk ro- 
tation into the system can be evaluated as 



nit) 



B^dRdz 



d^A$(*,t) . (20) 



We find that the ratio of J^^ to the initial total 

poloidal flux, J^z = Jo ° dR27rRBz = 27r, increases lin- 
early with time as expected by equation (^), reach- 
ing a maximum of ~ 4% when t = 4. So, the total 
injected toroidal fiux is still a small fraction of the to- 
tal initial poloidal flux. This information is useful for 
non-axisymmetric stability considerations for future 
studies. 

Figure |^ shows the spatial distribution of magnetic 
energy density and how it evolves with added twist. 
The quantity log^Q (i?^/Pc + 1) (i-e., a scaled mag- 
netic energy density) is shown for t = (top left) and 
4 (top right), where B'^ = B^ + B'^ and Pc = O.l. It 
is clear that essentially all the magnetic energy is en- 
closed by the plasma pressure. The expansion results 
in a clear collimation of magnetic energy along the 
z— axis. The two bottom panels show, at i = 4, the 
poloidal logiQ{Bp/Pc + l) (lower left) and the toroidal 



logxo (-B?/Pc + 1) (lower right) components, respec- 
tively. Note the dramatic increase of magnetic energy 
(mostly poloidal component) along the z— axis. 

4.4. Limitation on Increase of Twist 

We flnd that it is not possible to increase the twist 
beyond what is shown in Figure H. When more twist 
is added, the system evolves from a configuration with 
all field lines being tied to the disk to a new topol- 
ogy with some poloidal field lines close on themselves 
instead of connecting to the disk (i.e., forming an iso- 
lated "island" of poloidal fluxes or "plasmoid"). At 
this point, our method used to solve equation ([l3|) 
becomes unstable. The plasmoid formation is asso- 
ciated with the field lines that are in the FF regime 
(\J/ > ij/^). The transition to the plasmoid formation 
is sudden, giving a sense of eruptive behavior of the 
solution. However, as discussed below, this numerical 
behavior is associated with instability of the method 
but not necessarily related to instability of the MHD 
equations. 

Such eruptive behavior has been the subject of in- 
tense research in the solar flare community (Aly 1984, 
1991, 1995; Sturrock et al. 1995). The formation of a 
plasmoid from a single arcade (Inhester et al. 1992) 
is a loss-of-equilibrium bifurcation related to tearing 
instability of the current sheet which forms at the cen- 
ter of the arcade when it is sheared strongly (Finn et 
al. 1993). In addition, it has been shown that plas- 
moid formation can occur directly as a consequence of 
linear instability when multiple arcades exist (Mikic 
et al. 1988; Biskamp et a., 1989; Finn et al. 1992). 

This behavior of the solutions with increasing twist 
can be traced to a mathematical nature of the ellip- 
tical equation we are solving (e.g., equation (11) and 
(p^)). We can write this equation in a pseudo-time- 
dependent form, e.g., 

r)^ 1 dM^ 

— = ^-m + ----AR^^eM-{'^/'^c?]. (21) 

where a steady state is achieved progressively when 
d-^/dt ->■ 0. Here, A = SttPc/'^I is a positive 
constant. The computational domain can then be 
divided into three parts depending on the ratio of 
^/^c- Region I: large R. The pressure term 
^^—B?\dP/d^\^ dominates. Since it is negative, it 
drives ^ to exponentially for small ^ (i.e., d^/dt oc 
—Ci^, with Ci > 0). The solution in this region is 
quite stable. Region II: small R and small z (close to 
the disk), i.e., the FF region with a negligible pres- 
sure term. Region HI: small R but large z (along the 
rotation axis), where all three terms contribute. The 
eruptive behavior we found could occur in both re- 
gions II and HI, whenever the poloidal current term 



dH^ /d}^ exceeds a certain critical value. With a large 
current, magnetic fields tend to expand enormously 
and fill up the whole computational domain. 

An important question is whether this eruptive be- 
havior actually occurs in a system described with the 
full set of dynamical equations. Recent axisymmetric 
simulations using the full set of MHD equations by 
Ustyugova et al. (2000) and Goodson et al. (1999) 
indicate that field lines can reconnect and become 
open. Clearly, the rate of reconnection is enhanced 
by the artificially large resistivity in such codes; the 
exact role of resistivity in allowing reconnection when 
it might not otherwise occur is not completely clear. 
In Ustyugova et al. (2000), two regimes have been 
found after ~ tens of rotation periods of the inner 
disk: a hydromagnetic outflow from the outer part of 
the disk, and a Poynting outflow, which has negligi- 
ble mass flux but is dominated by the electromagnetic 
field along the rotation axis. 



5. DISCUSSIONS AND CONCLUSIONS 

We have shown that a static plasma pressure is fun- 
damental in shaping the overall magnetic equilibrium, 
despite the fact that the plasma pressure is exceed- 
ingly small (the maximum magnetic pressure over the 
plasma pressure is 10^^ in this study). This is essen- 
tially different from the pure force- free models (e.g., 
LB94; see however Lynden-Bell 1996). This differ- 
ence comes from the fact that the field lines that are 
being twisted the most expand the furthest so that 
they are most affected by the plasma pressure. This 
is opposite to the behavior in the solar flare models 
where most of the shear is concentrated around the 
0-point. The least amount of shear is applied near 
the 0-point in the accretion disk case. In regions 
where the plasma pressure is negligible, the physics 
of expansion is dominated by the force-free condition 
(J X B = 0). Their expansion is predominantly to- 
wards larger radii along an angle of ~ 60° from 
the z— axis (LB94). But for field lines affected by the 



plasma pressure, their radial expansion (away from 
the 2— axis) is slowed or stopped by this pressure. 
The twisting of those field lines causes them to ex- 
pand much more strongly along the z— axis due to 
the build-up of the toroidal field component B^. The 
build-up of Bfj, is non-uniform along those field lines, 
with more twist distributed on the ascending portion 
of the field lines. This is, again, a direct result of 
external plasma pressure. 

We have obtained magnetostatic equilibria up to a 
maximum of 4 turns of the inner most field line. For 
larger twist the solutions exhibit a change of topol- 
ogy and our method breaks down. The larger values 
of twist should be treatable by including the inertia 
term pdv/dt in the full set of dynamic MHD equation. 

An important question is the stability/instability of 
these axisymmetric equilibria. Instabilities in 3D in- 
volving flux conversion between toroidal and poloidal 
components is probably unavoidable and this effect 
will be important astrophysically: this is because the 
initial poloidal flux on the disk is probably too small 
to be responsible for the observed total magnetic flux 
in many astrophysical systems (Colgate &: Li 2000). 
The flux multiplication by the disk rotation/twisting, 
and subsequent flux conversion, could however gen- 
erate enough flux. 
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Fig. 1. — The distribution of poloidal fields on the disk between < R < Ront- Its ^ value is shown in the 
top panel, which increases from to 1 at the 0-point Rq = 0.01 and decreases back to at Rout = 0.02. Its 
corresponding \Bz{R,0)\ is shown in the middle panel, where it becomes zero as well as reverses its sign at Rq. 
The twist profile (oc R~^''^ is shown in the bottom panel. It is nearly Keplerian over large range of radii, but 
is flattened for i? — > and is approaching zero near the 0-point because the separation between the footpoints 
of the arcade is becoming zero. 
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Fig. 2. — The overall physical picture with two different regimes, separated hy ^ = ^c (the thick dashed 
line). The region labeled "plasma-pressure" (PP) is dominated by the plasma pressure with ^c ^ ^ ^ (the 
boundary at r = 1). A tiny fraction (=^c) of the total poloidal flux is in this region. The region labeled 
"force- free" (FF) is dominated by the magnetic pressure with ^c ^ ^ ^ 1 (the 0-point). Most of the poloidal 
flux (1 — ^c) is contained in this region. The key question we are addressing is the response of the ^ = ^c 
boundary to the twisting of field lines by the Keplerian disk. 
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Fig. 3. — The "evolution" of the poloidal field lines with increasing twist as solutions of equation (13) in a 
In r— spherical coordinate for t = 0,1,2 and 4, respectively. (We present the results in a smaller R — z plane for 
clarity.) The contours are displayed evenly in logarithmically spaced intervals (10"'^ < ^ < 1)- The outer field 
lines have expanded more strongly along the z— axis (from z ~ 0.1 to 0.3) than in the radial direction (from 
R^O.l to 0.15). 
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Fig. 4. — The poloidal current H{'^) (upper panel) and the twist A<I>(^) (lower panel) as a function of ^ for 
t = 1,2,4, (bottom to top curves), respectively, for the run shown in Figure y. The function H{^) is oc ^ for 
small ^ then flattens as ^ increases. The twist agrees with the input Keplerian profile nearly perfectly with 
indistinguishable differences. 
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Fig. 5. — The "movement" of three particular poloidal field lines, ^ = 10^^ (left), 10^'^ (middle), and 10~^ 
(right), as their twist increases from t = (solid), 2 (dotted) to 4 (dashed), using the results from Figure ^. 
The plasma pressure is important for ^ < ^c = 10~^. For field lines with \1' > ^c) they are force- free and 
expand along an angle 9 ~ 60° radially (the left panel). For field lines with ^ < ^ci their radial expansion is 
much "slower" than their vertical expansion. This results in the collimation around the z— axis. Again, these 
solutions were obtained in a In r— spherical coordinate with rmax = 1 but we present the results in a smaller 
R — z plane for clarity. 
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Fig. 6. — The distribution of cumulative twist along two field lines as a function of height z. The slope of 
the curves \d<\>ldz\ gives the pitch of the field line (twist per unit height). For \I' = 10~^ (dashed), its pitch is 
equally distributed between ascending and descending portions since it is force- free. For ^ = 3 x 10~^ (solid), 
its twist distribution is not uniform, with most of its total twist being distributed on the ascending portion of 
the field line. 



Fig. 7. — The magnetic pressure distribution at t = (top left) and t = 4 (top right). The quantity plotted is 
log^o(-^^/-fc + 1), where B^ = B^ + B'^ is the total magnetic pressure, Bp = B'j^ + Bl is the poloidal component 
and Pc = 0.1. The lower left and right panels are for logi^^Bp/Pc + 1) and logxo(^|/^c + 1), respectively, at 
t = 4. The dark blue color represents the region with log^Q(i?^ /i-*c + 1) ~ 0, i.e., B^ <^ P^. The added twist 
causes the fields to expand along the z— axis, with increased poloidal and toroidal magnetic pressures. Again, 
these solutions were obtained in a In r— spherical coordinate with rmax = 1 but we present the results in a 
smaller R — z plane for clarity. 
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